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Package-X: A Mathematica package for the analytic calculation of one-loop integrals 
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Package-X, a Mathematica package for the analytic computation of one-loop integrals dimension- 
ally regulated near 4 spacetime dimensions is described. Package-X computes arbitrarily high rank 
tensor integrals with up to three propagators, and gives compact expressions of UV divergent, IR 
divergent, and finite parts for any kinematic configuration involving real-valued external invariants 
and internal masses. Output expressions can be readily evaluated numerically and manipulated 
symbolically with built-in Mathematica functions. Emphasis is on evaluation speed, on readabil¬ 
ity of results, and especially on user-friendliness. Also included is a routine to compute traces of 
products of Dirac matrices, and a collection of projectors to facilitate the computation of fermion 
form factors at one-loop. The package is intended to be used both as a research tool and as an 
educational tool. 


Program summary 

Program title: Package-X 

Program obtainable from: CPC Program Library, Queen’s LTniversity, Belfast, N. Ireland, or 
http://packagex.hepf orge.org 

Licensing provisions: Standard CPC license, http://cpc.cs.qub.ac.uk/licence/licence.html 
Programming language: Mathematica (Wolfram Languange) 

Operating systems: Windows, Mac OS X, Linux (or any system supporting Mathematica 8.0 or 
higher) 

RAM required for execution: 10 MB, depending on size of computation 
Vectorised/parallelized? : No 

Nature of problem: Analytic calculation of one-loop integrals in relativistic quantum field theory 
for arbitrarily high-rank tensor integrals and any kinematic configuration of real-valued external 
invariants and internal masses. 

Solution method: Passarino-Veltman reduction formula, Denner-Dittmaier reduction formulae, and 
two new reduction algorithms described in the manuscript. 

Restrictions: One-loop integrals are limited to those involving no more than three propagator fac¬ 
tors. 

Unusual features: Includes rudimentary routines for tensor algebraic operations and for performing 
traces over Dirac gamma matrices. 

Running Time: 5ms to 10s for integrals typically occurring in practical computations; longer for 
higher rank tensor integrals. 


I. INTRODUCTION 

Many packages are available to assist with the evalu¬ 
ation of one-loop integrals that appear in higher order 
calculations of perturbative quantum field theory. The 
most widely used ones are the Mathematica packages 
FeynCalc[I], FormCalc[1| and the Fortran program 
Golem95[3]. These packages compute one-loop inte¬ 
grals using the Passarino-Veltman reduction algorithm [4] 
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(FeynCalc and FormCalc also feature a collection of 
routines designed to streamline the numerical computa¬ 
tion of a differential cross section; as such, they do sub¬ 
stantially more than to simply compute one-loop inte¬ 
grals). 

Nevertheless FeynCalc falls short in that it gives re¬ 
sults of one-loop computations in terms of basis scalar 
functions which cannot be evaluated on their own. In¬ 
stead, it is up to the user to supply their analyti¬ 
cal forms from an external source, or to link them to 
yet another package (such as FF[5], LoopTools[2], or 
OneLOop[S]). 

Moreover, one-loop integrals have many more applica- 
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tions than to calculate cross sections and decay rates. Ex¬ 
amples are the computation of ultraviolet counterterms, 
pole positions, residues, Peskin-Takeuchi oblique param¬ 
eters, electromagnetic moments, etc. Many of these ap¬ 
plications require the calculation of Feynman integrals 
at singular kinematic points such as at physical thresh¬ 
olds or zero external momenta. Since the Passarino- 
Veltman reduction algorithm typically breaks down at 
these points, it is nearly impossible to obtain results with 


Package-J calculates dimensionally regulated (d = 4 — 
( Pl ,..., PN - TOo , mi ..., m N ) = p 2e j jp- 


FeynCalc or FormCalc (Golem95 can give numeri¬ 
cal results). But, it is also at these points where compact 
analytic expressions exist. 

Although smaller-scale packages are available that are 
designed around a particular application (such as lool[7] 
and ANT 0). there is no general-purpose software that 
gives analytic results to one-loop integrals for all kine¬ 
matic configurations. In this regard, Package-J. serves to 
fill this gap. 

z) rank-P one-loop tensor integrals of the form 
1 •• • k^ p 

ml+ie\[{k+pi) 2 -rn\ +ie\ ■ ■ ■ [(k+pN) 2 -m 2 N -\-ie\ ’ 


with up to N = 3 denominator factors, and finds compact 
analytic expressions for arbitrary configurations of exter¬ 
nal momenta pi and real-valued internal masses irii. The 
functional paradigm of the Wolfram Language together 
with the supplementary trace-taking routines included 
in Package-J allows one to compute an entire one-loop 
diagram at once. All output is ready for numerical eval¬ 
uation and symbolic manipulation with Mathematica's 
internal functions. 

This article details the technical aspects of Package- 
X, and assumes familiarity in the use of the package. 
The application files are found at the Hepforge web¬ 
page http://packagex.hepforge.org, where the soft¬ 
ware will be maintained and periodically updated. In¬ 
cluded among the package files is a tutorial that provides 
an introduction, and a complete set of documentation 
files that becomes embedded with the Wolfram Docu¬ 
mentation Center upon installation which provides de¬ 
tails and examples of all functions defined in Package-J. 


II. STRUCTURE AND DESIGN OF PACKAGE 

The subroutines in this package belong to one of 
three Mathematica contexts organized as in Fig. |T] 
The module IndexAlg 1 contains the rudimentary tensor- 
algebraic routines and serves as the backbone of Package- 
X. OneLoop 1 contains the algorithms and look-up tables 
for the computation of one-loop integrals, and Spur 1 in¬ 
cludes the algorithms to perform traces over products of 
Dirac matrices and contains a catalog of fermion form 
factor projectors. 

The basic Package-J workflow for the computation of 
a one-loop integral consists of three steps: 

1. Call Looplntegrate to carry out the covariant ten¬ 
sor decomposition (section [TTT] ). 

2. Apply on-shell conditions and other kinematic 
relations with Mathematica' s built-in functions 
ReplaceAll (/.) and Rules (—>•). 



FIG. 1. Organization of functions into contexts as defined in 
Package-J. 


3. Call LoopRefine to convert coefficient functions 
into explicit expressions (section |Tv|). 

The reasoning behind the three-step design is as follows: 
kinematic configurations of external invariants Pi-Pj and 
internal masses to,; relevant to many physical applications 
occur at singular points of one-loop integrals, such that 
if they were applied after obtaining the general expres¬ 
sions, errors like 0/0 or 0 x ln(0) would inevitably occur. 
To avoid such errors and to facilitate the generation of 
compact results, LoopRefine uses algorithms depending 
critically on the kinematic configuration supplied by the 
user beforehand. 

Two other supplementary functions are provided to 
streamline computations involving fermions: 


Spur ( German for ‘trace ’) computes traces of Dirac 
matrices that may appear in the numerators of one- 
loop integrals (section VIII. 


Projector is used to project iernnon sen-energy 
and vertex form factors out of the loop integrals 
(section VIII). 


The algorithms used by these functions are detailed in 
the aforementioned sections below. 
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III. LOOPINTEGRATE: COVARIANT TENSOR DECOMPOSITION 


The evaluation of an integral is initiated with Looplntegrate, which carries out its covariant tensor decomposition 
in terms of scalar coefficient functions. For example (omitting the +ie), 


, 167T 2 


LoopIntegrate[k M k„k p , k, pi, mO, ml] : 

>-i 2e f d d k k^k u k p 

^ J (2n) d [k 2 - TOq][(/c + pi) 2 - m 2 ] 


{[pi][g]r vp B om + {[ Pl f}^ p B 


in > 


( 2 ) 


( 


l 

16tt 2 



Looplntegratefkpkjy, k, pi, p2, mO, ml, m2] : 

1 d d k k p k v 

(27r) d [k 2 — ml][(k + pi) 2 — to 2 ][(/c + P 2) 2 — m|] * 

MV'Coo + {\pi] 2 r v C n + {IpMV'Cu + {{P2] 2 T V C 22 


( 3 ) 


Here, I?ooi, S 111 , Coo etc. are coefficient functions that 
depend only on Lorentz invariants, Pi-Pj and m,;. Note 
that as indicated in the left hand sides, an overall con¬ 
stant is factored out of the natural integration 

measure p 2e to simplify the output. Each coeffi¬ 
cient function multiplies a totally symmetric tensor, de¬ 
noted {.. .} M "' in the notation of [9], containing products 
of external momentum four-vectors p/ and the metric 
tensor g pu . These tensors are generated by a Package -X 
internal function (inside IndexAlg 1 ), which utilizes Math¬ 
ematical built-in function Permutations. The time to 
generate the corresponding symmetric tensors grows fac- 
torially with the rank of tensor integrals. 

For integrals with high powers of contracted loop mo¬ 
menta, such as 


2e f d d k k a kP(k.k) 5 

^ J (2tt ) d [k 2 - ml\[(k+pi) 2 - mf\ ’ 

it is necessary to obtain explicit expressions of 
self-contracted symmetrized high-rank tensors like 
{bi] 6 M 3 } a/3M ' il ' 1/ppCT<TAA - K would be wasteful to first gen¬ 
erate the totally symmetric high-rank tensors, only to 
subsequently contract indices down to lower-rank sym¬ 
metric tensors. Instead, the contraction formulae 


(p fc u{bi] ni -"bivr N b] r r- pp 

N 

= -mm bir ••• [rf N brr- pp 

e =1 

+ (n* + i)(bfc]biP • • • biv]" N b]^ 1 }" 2 -^ (5) 

9 ^A[pir---\PN] nN [g] r Y^ p 

N 

= J2 Pi -Piibilbilbi ]" 1 ••• biv]" N b] r } P3 - pp 

i>3 

N 

+p 0 (d+P-2+^ n fc )(bi] ni •' • [PN} nN [gY- 1 }^ ^ , 
k 

( 6 ) 


are employed to carry out the self-contractions symboli¬ 
cally before converting any remaining symmetric tensors 
with free indices into explicit forms in terms of p/ and 
g 111 '. The time to construct self-contracted tensors in this 
way is reduced to follow a power law. 


IV. LOOPREFINE: REDUCTION TO ELEMENTARY 
FUNCTIONS 


Once the covariant decomposition is made, and any on- 
shell or kinematic conditions are applied, the final step is 
to feed the results into LoopRefine, which replaces the 
coefficient functions with explicit expressions. The basic 
algorithm followed by LoopRefine is as follows: 

Step 1: For each coefficient function (pvA, pvB, pvb, or 
pvC) encountered by LoopRefine, symbols for in¬ 
ternal masses are recorded (for Step 4), and the 
appropriate reduction routine (see corresponding 
subsections below) is called. 

Step 2: The reduction of C functions for more general 
kinematic configurations end with the scalar func¬ 
tion Co- If the C'o function is IR-divergent or has 
an explicit form that is sufficiently compact (as 
controlled by the option ExplicitCO), the explicit 
form is substituted. 


Step 3: All instances of the spacetime dimension d is re¬ 
placed by 4 — 2e, and Mathematical built-in func¬ 
tion Series is called to keep the leading terms in 
the e expansion. UV-divergences appear as 1/e 
poles, and IR-divergences appear as 1/e and/or 
1/e 2 poles. 


Step 4: Combine and simplify logarithms, organize the 
expression by the logarithms, and group the ’t 
Hooft parameter /^-dependent logarithm with the 
1/e pole in the expression (see section VII. 


In the following subsections, the algorithms and ac¬ 
companying formulae used by LoopRefine to reduce the 
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coefficient functions are summarized. It should be noted 
that nearly all algorithms are drawn from the 2005 pa¬ 
per by Denner and Dittmaier I5|, and will be referenced 
henceforth as [DD]. The only formulae not taken directly 
from their paper are those for the auxiliary b& functions 
in Section |IVB (which is only a slight modification of 
the reduction formulae for B functions), and those of 
two additional algorithms for the reduction of C func¬ 
tions in special kinematic configurations (Cases 2 and 4 
in section IV C). 


A. Reduction of A and B functions 


where 1/e = 1/e—7B+ln(47r), and H n is the ?r th harmonic 
number. 

The f?o...oi...i functions, with at least one pair of 00 
indices are obtained iteratively in terms of those with 
fewer number of 00 indices using (eqn 4.5 of [DD]): 

= 2 (n + l) [(~l) n+1 A°^°,( mi ) 

2r n ^ 2(r—1) 

+ (p 2 - to 2 + ml)B 0 „,o 1...1 (p 2 ;m 0 ,mi) 

2(r—1) n + 1 

+ 2p 2 B Q ,„ Q 1...1 (p 2 ;to 0 ,toi)] , r> 1 (8) 

2 ( 7 — 1 ) 71 + 2 


The Passarino-Veltman coefficient A functions are sim¬ 
ple enough to be obtained by direct integration (eqn 3.4 
of [DD]): 


^ (m °) - nr + 1)! (l + ln( ^ } + Hr+1 ) ’ (7) 


2 r (r + 1) 


Then the integrals (with no 00 index pairs) are 

obtained by explicit integration over the single Feynman 
parameter in ( |B2[ ). Results are given in (eqn 4.8 of [DD]), 
but the form that tends to generate most compact expres¬ 
sions is 


S 1...1 (p 2 ; mo, mi) = 


(-l) r 


n + 1 Le 

L^J 

E 

k =0 

LfJ 

+ E 

k =0 


1 LL 2 n 2 OP 

Eln(^)+V^-V 

/= ' mr ' ^ n .4- 1 ' 


n + l 

k =0 j—0 


n — k\ f p 2 + to g — mf 

j ) 1, v 


( A(p 2 ,TOg,TO?) 

1 4 (p 2 ) 2 


n + 1^ ^p 2 + TOq to 2 ^ ” +1 2k ( A(p 2 , TOq, to 2 ) ^ k ^ ^ TOq 


2p 2 


\ 4(p 2 )2 


n + 1 A /p 2 + TOq - ?n 2 A " / A(p 2 , TOq, mf) 

2fc +1/ \ V 7 V 4(p 2 ) 2 


A(p 2 ; TOq, TOi) 


(9) 


Here A (a, 6, c) = a 2 + & 2 + c 2 is the Kallen function, implemented as KallenA[a, b, c], and A(p 2 ;?no,mi) is the 
abbreviation 


A(p 2 ; to 0 , TOi) = 




TOq,TOJ 


■ In 


2toqTOi 


—P 2 + TOq + TO 2 — \/A(p 2 , TOq, TOi) 


■ IE 


( 10 ) 


implemented as DiscB[s, mO, ml] . In order to access (p 2 ; too, toi) at its singular points, a limiting procedure would 
need to be made at runtime in order to avoid errors such as 0/0 or 0 x ln(0). While Mathematical s function Limit 
can eventually generate an expression, computation time is long, and output expressions are always unwieldy. Instead, 
a catalog of explicit expressions (also obtained by direct integration) of H 1...1 at all its singular points (see Table [i]) is 
included in the source code. They may be accessed directly within Package-X using LoopRef ine[pvB[0, n, s, TOq, toi]]. 


B. Reduction of auxiliary // functions Veltman f/ functions |10| . For example, 


In covariant gauges, the propagator for massless vector 
fields 


iD^(k) 



( 1-0 




( 11 ) 


i \ —1 2e f d d k k^khk 9 

M J (2 nY [fc 2 ] 2 [(fc+p) 2 - to 2 ] 

= {\p][9}r Vp bU+{\p} 3 V VP biu- 


contains a gauge term that leads to an additional fac¬ 
tor in the denominator of one-loop integrals. Package-X 
can handle such propagators inside bubble integrals, with 
the coefficient functions given by the auxiliary Passarino- 


The reduction formulae for these functions essentially 
mirror those for the standard B functions. Auxiliary 
6 q 01 1 functions with at least one pair of 00 indices are 
iteratively determined in terms of functions with fewer 
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00 index pairs using 


^o...o 1...1 (P i m ) — 


-1 



[-B0...0 (p 2 ; 0 ,m) 


+ (P -m )b%_ 0 1...1 (P ;"i) 


+ 2 P 2fo o...o 1...1 (p 2 ; TO )] > r>l, (12) 

2(r-l) n + 2 

and the 6^ 1 functions with no 00 index pairs are ob¬ 
tained by direct integration over the single Feynman pa¬ 
rameter in (B3). The integral is finite if n > 1, with the 
result 


v ^,( p2;m ) = 

(- 1 ) 


n— 1 r 1 ™^ 


-- + E — 

n z ' tj — , 


n z —' n — k p 2 — m 2 

fe=i 


m 


p 2 — m 2 


p 2 — m 2 


In 


to 2 — p 2 


2 2 \ ^ 
p — m 


■ 


If n = 0 (a case that is not met in practice since the gauge 
part of the spin-1 propagator guarantees two powers of 
momenta in the numerator), the auxiliary b& function is 
IR-divergent. 

Explicit expressions at the various singular points of 
b{ A (see Table [l| are included in the Package-X source 
code. 


C. Reduction of C functions 


In the reduction formulae below, the following kine¬ 
matic abbreviations are used (which differ slightly from 
[DD] by numeric factors): 


fj =p 2 -m 2 + ml, j = { 1,2} 


z= ( Pi Pip) 
\P2-Pl P 2 J 

q 2 =p\+p\- 2pi.p 2 

det Z = \\{q 2 , p \, pi) 


(Gramian matrix) 


(13) 


f pI ~p\p?\ 

\-pl.p 2 p\ ) 


_ ( P 2/1 -Pl-P 2 f 2 \ 
0j \-Pl-P2fl +P 2 lh) 


(cofactor matrix) 


3 = (1, 2} 


Furthermore, hatted indices on coefficient functions (e.g. 
B~ k0 01 1 ) indicate the removal of those indices. Co¬ 
efficient B functions derived by canceling denominators 
from three-point integrals are abbreviated by 

B...(Di) =B..Xpl;m 0 ,m 2 ) (14) 

B..XD 2 ) =B...(pl;m 0 ,m 1 ). (15) 


If the denominator (fc 2 — Wq) _ 1 independent of an exter¬ 
nal momentum vector is cancelled, a shifted form of the 
B function is used: 



(16) 


The reduction of coefficient C functions is significantly 
complicated by its numerous singular points. Although 
the standard Passarino-Veltman reduction algorithm is 
applicable at almost all points (Case 1 below), differ¬ 
ent formulae are needed to handle the various singular 
cases (Cases 2-6). LoopRefine identifies the nature of 
the kinematic configuration and applies the appropriate 
reduction method. 

Cases 1, 3, 5 and 6 are taken from [DD], Note that 
since the emphasis of [DD] is on numerical stability and 
not on generating analytic expressions, the algorithms 
presented there do not automatically give compact ex¬ 
pressions. The algorithm under Case 2 is new, and while 
technically it is covered by Case 1, it leads to more com¬ 
pact expressions. Furthermore, an algorithm to handle 
the reduction at physical thresholds (applied in Case 3 
below) is not completely covered by [DD]. This gap is 
filled by the formulae under Case 4- 

The arguments of the coefficient C functions are or¬ 
dered differently in Package-X as compared to those used 
by other authors. See Appendix [A] for details. 


Whenever n\ > n 2 the invariance property 
B 0...0 1...1 2...2 (Do) = B 0...0 1...1 2...2 (Dq) 


(17) 


is used to keep the number of terms in the sum to a 
minimum. Cases 2 and 4 require expressions for the B 
functions with the number of 00 index pairs continued to 
r = — 1. Details of this function are found in Appendix 

m 

Finally, formulae for Cases 1 , 3 and 5 below contain 
explicit dependence on spacetime dimension d = 4 — 2e 
appearing in denominators of certain prefactors. In the 
course of reduction, the O(e) part multiplying any lower 
coefficient functions combines with their UV 1/e pole^J 
and gives rise to finite polynomials in kinematic variables. 
Although this can be automatically handled by Series 
at Step 3, the reduction algorithm performs much faster 
if these polynomials are explicitly supplied. They are 
obtained by integration over the Feynman parameters as 
described at the end of Appendix [B] 
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Case 1: det Z ^ 0 

At non-singular kinematic configurations with det Z ^ 0, the original [1] Passarino-Veltman reduction formula is 
used (eqns 5.10, 5.11 of [DD]): 


C 


1 


0...0 1...1 2...2 — 


2 det Z 


E Zj k 


k=l 


fin k , S jk B 0...0 1...1 ( Dk ) — B 0 0 1...1 2...2 { Do ) 

0...0 1...1 2...2 (A>) -( ^k 0...0 1...1 2 ... 2 J 7 ^ (18) 


C 0...0 — 


^ 2 ? .~\ [Bo...o (A>) + 2toqCo...o i + / 2 C 0...0 2 J , 


2r+2 ni-1 n 2 


2(d — 4 + 2r) 


r > 1 


1 k = 2 

where k = < .In the first equation, j = 1 is taken, although the choice j = 2 would give equivalent results. 

2 , k = 1 


If m = 0 with ri 2 > 0, then the relation (B5) is used and the first equation is applied. 


Case 2: Ellis-Zanderighi triangle 6 


Coefficient C functions for which arguments are (toq, s, m^; m- 2 ,0, too) —or an equivalent permutation thereof—are 
already covered by Case 1. However, final expressions obtained from it tend not to give the most compact formulae 
for this kinematic configuration. More compact formulae are obtained by directly integrating over the Feynman 
parameters in (B41; see Appendix [C] for derivation. It is of note that the corresponding scalar function Cq is the 
IR-divergent three-point function, ‘triangle 6’, as classified by Ellis and Zanderighi m- In Eqs. |19| ) and ( |20| , it is 
assumed that at least one of r, to or n ,2 is nonzero. 


C'.o...o..i...i,. 2...2 (ml, s, ml; m 2 ,0, m 0 ) = 

(-l) ni ni!(n 2 + 2r - 1)! 
2 (to + n 2 + 2r)! 


C 1...1 (toq, s, ml;m 2 , 0, to 0 ) = 


(id - 2e(7? ni _|_ n2 _|_2T- FT n2 _|_ 2r—l)) By., .o^ i ... i (s, too , to-2 ), 

71 2 

n 2 ^ 0 or r ^ 0 (19) 


(-l) ni Co(too, s, to|; to 2 , 0, too) - + e(ff^ 1 + H^)^B 0 ...o (s; m 0 , m 2 ) 


where Hj' 1 is the n th harmonic number of order r. If the arguments take the form (s, Toq, to|; 0, m|, Too), then the 
identity (B5) is applied, and the equations above are valid. 


A different formula is needed if the off-shell momentum s is in the third position: 


C 0 ...o 1...1 2. ..2 (to2,TOo,s;to 0 ,to 2 ,0) = 


(-l) r 


TO + n 2 


2 r 


k—0 


e r; (!+ 


2e 


ni + n 2 + 2 r 


I-B.o.. .o..i.. .i,(s; m 0 , m 2 ) (20) 

il +fe 


To apply Eqs. (19) and (20) above, explicit forms of the scalar B function with the number of 00 index pairs taken 
to r = — 1 is occasionally needed. These functions are discussed in Appendix [Ej 


1 For the argument that they are only of UV origin (and not IR), 
see the argument in Sec. 5.8 of [DD] 
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Case 3: det Z = 0 but Xoj ^ 0 

With det Z = 0, the primary reduction formulae are rearranged to give: (eqns 5.38 and 5.40 of [DD]) 


C’o.-.o — 


1 


d + 2r — 3 

2 r 2r — 2 

2 

X ^ ' Z n j (1 Smj )B Q...0 1 (Dm) 


Bq...o (Do) — ffl'oCo...o ) + 


1 

2 (d + 2 r — 2>)Zi~i 


S (■ 


^km&nl d kjSr, 


- B 


j =i 


i.cr-o 

2r —2 


(A>) 




— B0...0 (D n ) + B o...o (-Do) + fnCo...o J j r > 0 


1 z 

= E] -^jfe — B a^o /v i_i / ^_2 / (£ ) o) — 2 rikC l 


fcOO 0...0 1...1 2...2 


(2i) 

The value of j chosen (1 or 2) is the one for which the corresponding Xoj is non-vanishing. If both elements are 
vanishing, then Case 4 is applied. Note that the second relation is valid even when either rq = 0 or n 2 = 0. In 
particular, when r = n\ = n 2 = 0 the final term in that relation vanishes, and leads to the reduction of the scalar Cq 
function in terms of scalar Bq functions. 


Case 4'- vanishing det.Z and Xoj 

When the physical threshold (corresponding to Xoj = 0 for both j = {1,2}) coincides with the boundary of 
the physical region (det Z = 0), then Cases 1—3 are inapplicable. For this kinematic configuration, the reduction 
formulae in [DD] eqns (5.49) and (5.53) can be used provided at least one element of 

4 m 2 0 p 2 2 - /f -2 ml(p\ + p 2 2 - q 2 ) + fif 2 

-2 ml(p\ + v\ - q) + / 1/2 4 mlp\ - ft 

is non-vanishing. However, no reduction methods are presented in [DD] that are valid when all four elements of X. t j 
are vanishing, because an expansion around that point is not knowrrl This exceptional configuration is needed for the 
computation of elastic form factors at zero momentum such as electron g — 2. To fill this gap, a new set of reduction 
formulae are used that is valid regardless of the form of X,j , provided at least one of p 2 , p\ or q 2 is non-vanishing. 
These formulae are derived in Appendix [D] 

If pi ± 0, 



^ _ (_i)ni+n2 ^ 

^.0...0„1...1„2...2. — -^ 

3 =0 

(-ir 


a 


n 2 -J 


»l!(^2 - j)'- 
(ni + n 2 - j + 1)! 


E 

k—0 


(-ay- k (-l) k B 0 ... 0 (Di) 


2i— 2 k 


+E — 


k—0 


n 2 — j + k + 1 \ k 


n\ 


(1 - ay +1 (-l) n2 B^3 na (4) - (-a) j+1 B^^i (4) \, (22) 

2r-2 n 2 +k +1 2r-2 n 2 + k +1 


where a = —q 2 + p 2 +p|/(2p|). 

If P 2 = 0, then det Z — 0 implies q 2 = p\ , and the formula 


_ (-l) ni+1 y- / ni \ 

f O' 0 ' 1 2 2. 2( „ 2 . n 2 ^l ) /; 

n 1 n 2 


k—0 


0...0 1...1 ( d 2 ) 

2, -2 n 2 + k + l 


(23) 


is used. If p 2 = p 2 = q 2 = 0, then these formulae are inapplicable and Case 5 is needed. Note that when r = 0 in 
either (221 or (23), the B functions continued to r = — 1 are needed (see Appendix [E}. 


2 A. Denner, private correspondence 


















Case 5: All elements of Z vanishing 


If all external invariants are vanishing p\ = p 2 = q 2 * 
[DD]): 


0, then the following are applied (eqns 5.62 and 5.63 of 



d + 2(ni + tt -2 + r — 1) 

1 r 


Dq...O 1...1 2...2 (Do) + ITIqC 



r > 1 


— ~r S n k oBi...i {Dk) — B i .1 2 ...2 (Do) — 2rikCj 


kOO 1...1 2...2 


(24) 


k 71 1 71 2 n k n l 71 2 ni n-2 

In the second equation, the index k is chosen such that fk is non-vanishing. As in Case 3, the second relation is valid 
for vanishing n\ or 712 , and is used to reduce the scalar Co function to Bq functions for n\ = 712 = 0. 


Case 6: All elements of Z and fk vanishing 


Finally, if also the fk s are vanishing, the following formulae are used (eqns 5.71 and 5.72 of [DD]): 



Bk 0...0 1...1 2...2 (Do), 

2r-2 »i n 2 


r > 1 


— —2 + 2 tii + 2 t 7 2 )Coo 1...1 2...2 — -B1...1 2...2 (Dq) 


(25) 


Equivalent results are obtained for k = 1 or 2 in the first equation. The choice k = 1 is used in Package -X. Note that 
when r = Tii = 712 = 0, these relations together permit the reduction of the scalar Co function in terms of scalar Bq 
functions. 


V. THE SCALAR C 0 FUNCTION: ANALYTIC 
EXPRESSIONS AND NUMERICAL 
IMPLEMENTATION 


The algorithms for the reduction of coefficient C- 
functions for which det Z / 0 (Cases 1 and 2 in the 
previous section) end with the UV-finite scalar func¬ 
tion Cq(pi, p%,d 2 1 m 2 ,mi, mo). To complete the compu¬ 
tation of the one loop integral and to make the final out¬ 
put useable, the scalar function must be replaced. For 
this purpose, a complete catalog of analytic expressions 
for Co where det Z ^ 0—each one obtained by direct 
integration—is included in the source file (see Table i]). 
Many such expressions are scattered throughout the lit¬ 
erature. The general formula with non-zero kinematic 
variables appears in [42'. All IR-divergent three-point 
formulae are given in mi, and some special cases appear 
in unpublished notes m 

Although a complete catalog of analytic expressions 
of Co is available, not all cases are automatically sub¬ 
stituted by LoopRef ine at Step 2. The functions that 
are substituted are only those that are IR-divergent (to 
faithfully display the 1/e poles in the final output), and 
those for which a sufficiently simple/compact expression 
is known. For more complicated finite cases, LoopRef ine 


simply output^] pvC0[. ..], with the function itself imple¬ 
mented numerically, (summarized below). The reason for 
this design choice is as follows: 

Firstly, in cases for which no simple form is known, the 
general formula [12] in terms of 12 dilogarithms would 
have to be given. This expression for Co alone would 
occupy a very large part of the output overwhelming the 
remainder of the expression, thus defeating the original 
purpose of producing compact expressions. Secondly, the 
dilogarithm function is computationally very expensive. 
When numerics are required, a brute-force evaluation of 
all the dilogarithms is grossly inefficient, leading to ex¬ 
cessively slow numerical evaluations. 

The main features of the code for the rapid numeri¬ 
cal evaluation of the three-point scalar function for real 
masses and external momenta are as follows: 

• The imaginary part of Co in the physical region 
(defined by A(q 2 ,p\,p\) > 0) is obtained by ap¬ 
plying Cutkosky’s rule, and with a straightforward 
continuation into the unphysical region (defined by 
iPiiPi) < 0) [ISUS]. Its computation requires 
the evaluation of a single logarithm. 


3 If the explicit analytic form is desired, the option 

ExplicitCO —> All can be supplied to LoopRef ine. 
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B-functions — Section IV A 


Bi...i(0;0,0) 

B 1...1 (0; 0, mi) 
B 1...1 (0; mo, 0) 


Bi... i(p a ;0,0) 

-Bi...i(0; mo, mo) 
0 , mi) 


B-functions with r 

Si...:l(0;0,0) 

B 1...1 (0; 0, mi) 
-Bi...i(0;m 0 ,0) 


-1 


Appendix IEI 


Bi...i(p 2 ;0,0) 
Bi...i(0; m 0 , mo) 
Bi...i(mi; 0, mi) 


Auxiliary -functions — Section IV B 

W. .1(0,0) OW" 

fof !(0,m) b\ A {m 2 ,m) 


— Section V 


Scalar C-functions 

C 0 (0, 0,0; 0,0,0) 

C 0 (0,0,g 2 ; 0,0,0) 
Co(0,0,ml;m 2 ,0,0) 
Co(0,0,g 2 ;m 2 ,0,0) 
Co(0, 0, q 2 ; 0, 0, mo) 
C 0 (0,0, m 0 , m 0 , m 0 ) 


C 0 (0,0, q ;0,m o ,mo) 

Co(0, 0, g 2 ; 0, mi, m 0 ) 
C 0 (0, 0, q 2 ; mo, m 0 , m 0 ) 
C 0 (0,0,g 2 ;m 2 ,m 0 ,mo) 
C 0 (0,0,g 2 ;m 2 ,mi,m 0 ) 
Co(Pi, 0, g 2 ; 0,0, 0) 


Bi...i(mj; mo, 0) 

Bi...i(0;m 0 ,mi) 

-Bi...i(p 2 ;0,mi) 


-Bi...i(p 2 ;m o ,0) 

-Bi...i((m 0 +mi) 2 ; m 0 , mi) 
Bi...i((mo-mi) 2 ;mo, mi) 


Bi...i(mo; m 0 ,0) 
-Bi...i(0; m 0 , mi) 
Bi...i(p 2 ;0,mi) 


-Bi...i(p 2 ;m o ,0) 

Bi...i((mo+mi) 2 ;m 0 , mi) 
-Bi...i((m 0 -mi) 2 ; m 0 , mi) 


Co(mo, 0, g 2 ; 0, 0, m 0 ) 

Co(0, m 2 , g 2 ; m 2 , 0, 0) 

Co (m§, 0, ml; m 2 , 0, m 0 ) 
Co(pl,pl,q 2 ;m 2 ,mi,m 0 ) 
Co(0,ml,g 2 ;0,m 0 ,m 0 ) 
Co(0,p!,g 2 ;m 0 ,mo,mo) 


C o (0,pl,g 2 ;m 2 ,0,0) 

Co(p 2 , 0, g 2 ; m 2 , mi, m 0 ) 
C o (ml,ml,g 2 ;0,0,m 0 ) 
C o (m5,pl,ml;m o ,0, m 0 ) 
C o (m§,pl,ml;m 2 ,0, m 0 ) 
C 0 (p?,pl,g 2 ; 0,0,0) 


TABLE I. Special kinematic cases of the Passarino-Veltman coefficient functions B, b 5 and C for which explicit expressions are 
included in the source file OneLoop. m. Further information for these functions is found in the indicated sections 


• The real part of Cq requires evaluations of only 
the real part (in the physical region) or only the 
imaginary part (in the unphysical region) of the 
dilogarithm, but not both. Calling PolyLog would 
lead to needless computation of both parts by the 
Mathematica Kernel. Following m , the real and 
imaginary parts of the dilogarithm function are im¬ 
plemented separately. 

• For the real part of Cq in the physical region, the 
+is prescription is irrelevant (since it influences 
only the imaginary part which is anyway evaluated 
using Cutkosky’s rule). Then, either the arguments 
of the 12 dilogarithms come in complex-conjugate 
pairs (for which the real part of the dilogarithms 
are identical and are added reducing the number 
of dilog evaluations), or the arguments are purely 
real (for which the real parts of the dilogarithms 
are rapidly evaluated using real arithmetic). 

• The code is compiled to the Wolfram Virtual Ma¬ 
chine (using Compile), leading to a substantial 
boost in computation speed. 

In the physical region, up to a 200-fold increase in 
speed is achieved compared to brute-force Mathematica 
evaluation by the Kernel. In the unphysical region, up 
to a 20-fold increase in speed is obtained. If the op¬ 
tion Compilation!arget —> “C” to Compile is set, its per¬ 
formance rivals that of the Fortran implementation in 


LoopTools, with Package-X generating results approx¬ 
imately twice as fast. 


VI. HANDLING THE -Me PRESCRIPTION AND 
SIMPLIFYING LOGARITHMS 

The +ie prescription appearing in the denominators 
of propagator functions enforce causality in the time- 
ordered Green functions of a relativistic quantum field 
theory. In one-loop computations, it determines the 
branch on which the logarithms are to be evaluated. All 
output expressions of LoopRefine observe the +ie pre¬ 
scription and are consistent with the analytic conventions 
of the built-in Mathematica functions Log and PolyLog, 
which are 


Log[x] —i lim ln(x + ie ), and 

E-S-0 + 


PolyLog[2,x] —> lim Li 2 (x — is). 

£->0+ 


Because Package-X assumes real external invariants and 
internal masses, almost all analytic formulae can be ex¬ 
pressed compactly in terms of the built-in functions. 

Whenever LoopRefine generates a logarithm contain¬ 
ing the ratio of two internal masses, the ratio may be 
flipped to bring the logarithm to ‘canonical form’, e.g. 


Log 


‘ml 2 ' 
.mO 2 . 


-Log 


‘m0 2 ‘ 

.ml 2 . 


(26) 
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Since internal masses are assumed to be positive real, 
this is allowed, and helps to keep the logarithmic parts 
compact. 

In the course of reduction, regardless of whether the 
final expression is divergent or finite, multiple logarithms 
of ratios of several scales with the’t Hooft parameter p 2 
are typically generated, e.g. 


a Log 


‘ /iR 2 ' 
. —s . 


+ bLog 



+ c Log 



(27) 


It is found that by consistently keeping p 2 in the numera¬ 
tor, the +ie prescription is always observed - even when 
the other scales are external invariants that may become 
time-like. The p 2 from each logarithm are brought into a 
single logarithm by forming the ratio with a variable that 
is known to be positive (which were recorded at STEP 1): 



a Log 



+ (a + b + c) Log 



+ cLog 



(28) 


That way, the coefficient- (a+b+c) in this example—of 
the ^-dependent logarithm always matches that of the 
1 /e pole elsewhere in the expression, and are grouped 
before presenting the results. If the final expression were 
in fact finite without a 1/e pole, the coefficient would 
cancel exactly. 

In more complicated cases, expressions cannot be given 
compactly assuming a universal sign for the infinitesi¬ 
mal imaginary part. In this case, Mathematical s built-in 
functions Log and PolyLog are unsuitable. For this pur¬ 
pose, two new analytic functions are defined in Package-1. 
(within OneLoop 1 ): 


Ln[x,a] —> lim ln(cc + iae ), and 

e->0+ 

DiLog[x,a] —► lim Li 2 (a; + iae). 

e->0+ 


The (real part of the) second argument a controls the side 
of the branch on which these functions evaluate. A simple 
example that uses DiLog in its output can be found by 
running 


LoopRef ine[pvC0[0, m 2 , s, 0, m, m]] . 

VII. SPUR: COMPUTATION OF TRACES OF 
DIRAC MATRICES 

To assist in the evaluation of one-loop integrals with 
internal (closed) fermion lines, Package-1 includes the 
rudimentary function Spur (inside the module Spur 1 ) to 
evaluate traces over products of Dirac gamma matrices 
appearing in numerators. The function Projector helps 
to handle loop integrals with open fermion lines and is 
described in the next section. Because the primary func¬ 
tion of Package-1 is to compute loop integrals, with the 


computation of traces being a secondary feature, only 
a cursory description of the algorithms are given in the 
following two sections. 

As with the rest of the algorithms in Package-1, the 
calculation of traces is rule-based at its core, and bears 
some resemblance to that of a much earlier Mathematica 
package Tracer|T7]. However there are a number of 
differences listed below that lead to greater computation 
speed. 

• Throughout the evaluation of the trace, expres¬ 
sions can grow very large containing many terms. 
Groups of terms are temporarily enclosed within a 
List to prevent the Mathematica kernel from au¬ 
tomatically simplifying the large expression at each 
step of the computation process. 

• While more complicated trace formulae such as 
those for products of numerous gamma matrices are 
recursive, non-iterative rules are used for simpler 
tasks such as for collecting 75 and Pl/Pr within 
each term. 

• Products of gamma matrices with repeated Lorentz 
indices (such as 7 /i 7 1/ 7 p 7 /i ) are related to products 
with fewer gamma matrices. With more gamma 
matrices interposed between contracted matrices, 
the number of terms in the identity grows. On 
account of the cyclic property of the trace, these 
contraction identities may be applied in one of two 
directions. Additional rules are included so as to 
apply the identity in the direction with fewer num¬ 
ber of interposed gamma matrices. 

• Traces that multiply 75 are tagged differently to 
set it apart from those without it. This way, rules 
for computing traces with 75 and those without 75 
are separated, and saves some time when the kernel 
searches for the appropriate rules. 

When compared to the other Mathematica packages 
FeynCalc and Tracer, Package-1 generally gives re¬ 
sults around 10 times faster. As an example, the trace 

Tr [(# - - f 2 + m)Y(gLpL + gRpR){ft -f 2 + m) 

YigLpL + gRpR)(fi + m)"i >1 (g L P L + gRpR)\ (29) 

was calculated with each package and computation times 
were recorded (with Timing). The results on a 2.93 GHz 
Intel i 7 processor are: 

Package-1 0.096 s 

FeynCalc 8.2.0 1.03 s 

Tracer 1.1 0.81 s 

Part of the motivation for refining the trace-taking al¬ 
gorithms is due to the inclusion of fermion projectors 
described in the next section. When a Projector is in¬ 
cluded inside Spur, the number of terms within the trace 
is increased to an extent that a noticeable slowdown is ob¬ 
served. However, with the refinements described above, 
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projections onto form factors are nearly instantaneous on 
a modern computer. 


VIII. PROJECTOR: PROJECTION ONTO FERMION 
FORM FACTORS 


Package-X does not directly handle expressions involv¬ 
ing open fermion chains that are relevant for fermion self 
energy and form factor calculations. In order to provide 
some support for such computations, Package-X comes 
equipped with a set of projectors. The projectors permit 
the projection of a loop integral with an open fermion 
line onto specific form factors functions. 

For example, the one-loop expression for the off-shell 
fermion self-energy function takes the form 


1(f) = M 2e 


d d k M (k,p) 

(2n) d [k 2 — m 2 ] [(k + p) 2 ] 


(30) 


where M(fc,p) is a Dirac matrix structure that depends 
on the integration variable k and external momentum p. 
Parity conservation and Lorentz covariance allow I to be 
written in the form 


1(f) = A(p 2 )f + B(p 2 )m , (31) 

where the form factors A and B depend on Lorentz in¬ 
variants p 2 and to 2 only. By multiplying the appropriate 
projectors 

T^(p,m) = and T [B] (p,m) = ^ 

with the numerator of ( |30| , and taking the trace, the 
form factors are obtained: 

2 2e r d d k TV[M (k,p) P^ A \p,m)] 

' [P ’ P J (2n) d [k 2 - m 2 ][(k + p) 2 } 

p( 2 \ _ 2e /' d d k Tr[M(k,p)pW(p,m)] 

[P ’ P J (2n) d [k 2 - m 2 ][(k + p) 2 } ' 

The trace over the projectors convert the expressions into 
ordinary tensors integrals that are readily computed with 
Package-X. 

A large set of pre-programmed projectors for off- 
shell self energy functions and on-shell scalar- and 
vector-vertex functions in various bases (L/R-chiral or 
Vector/Axial-vector) are available (as Projector) to 
streamline the computation of such integrals. These pro¬ 
jectors are generalizations of those used in m for the cal¬ 
culation of lepton anomalous magnetic moments, and in 
m for dipole moments. A comprehensive list of available 
projectors is given in the built-in documentation files. 


Table [D The reduction routines for A and B coefficient 
functions and the basis functions were compared 

against another (unpublished) computer program devel¬ 
oped by Huaike Guo. The reduction of C functions for 
Cases 1, 3, 5 and 6 were checked against explicit for¬ 
mulae for the low rank functions listed in [DD]. Each 
Cq scalar function was derived by hand and compared 
against explicit formulae in the literature where they ex¬ 
ist 0 QH H3 ED]- In cases where they did not exist, 
the analytic expressions were checked by comparing with 
the results of numerically integrating the corresponding 
Feynman parameter representations given in Appendix 
0 

Finally, as a combined check of the various algorithms 
in Package-X , the following well known physical quan¬ 
tities were computed and verified: H —► gg and 77 
standard model decay rates [21 j, electron g — 2 , and the 
neutrino electric and magnetic moments [22) . Each was 
found to be in agreement with literature. 

There are a number of important limitations of Pack¬ 
age-X , listed below, that guides its current line of devel¬ 
opment. 

1. An analytic series expansion of the loop integral in 
kinematic variables is not generally possible. Cur¬ 
rently, the only available method is to use Math- 
ematica's Series on the output of LoopRefine. 
However, if the result of loop integral contains spe¬ 
cially defined function like pvCO, then Series will 
not work. Given that much information about a 
loop-integral can be gleaned from its expansion, the 
omission of this feature is most conspicuous. 

2. Package-X currently supports loop integrals with 
up to only three denominator factors. But, as the 
number of denominator factors increases, so does 
the complexity of their analytic forms. Thus, it 
would not be so practical to work with such ex¬ 
pressions for higher-point functions even if Pack¬ 
age-X were to provide them. However, at special 
kinematic points such as at zero external momenta 
or at thresholds compact expressions could be ob¬ 
tained. 

3. Gamma-5 is implemented naively in dimensional 
regularization. This means that the VVA or AAA 
three-point functions may not automatically satisfy 
Ward identities appropriate to the physical prob¬ 
lem. However, the versatility of Package-X makes 
it easy to apply Adler’s method [23] (see also El) 
to enforce the Ward identities. 


IX. CROSSCHECKS AND FURTHER 
DEVELOPMENT 


The verification of loop integrals obtained by Package- 
X is divided into two parts: checking the reduction algo¬ 


rithms in Section IV and checking the basis functions in 


4. Loop integrals with open fermion chains are not 
directly supported. As explained in Section |VIII[ 
there is no way to input an open string of Dirac 
matrices. Instead, the computation of fermion form 
factors can be done by projecting out the needed 
form factors. 
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Appendix A: Conventions 


For reference, the conventions for spacetime quanti¬ 
ties are summarized in Table El Conventions for the 
Passarino-Veltman functions are displayed in Table m 
Note that a sliglitly-unconventional ordering and form 
for the arguments of the Passarino-Veltman functions is 
taken. However, this choice makes the invariance prop¬ 
erty under their pairwise interchange clear: 


Co{pl,pl,q 2 ,m 2 ,m 1 ,m 0 ) 

= C 0 (pl,pl,q 2 ,m 1 ,m 2 ,m 0 ) 

= C 0 {q 2 ,pl,Pi,m 0 ,m 1 ,m 2 ). (Al) 


Appendix B: Feynman parameter integral 
representations of Passarino-Veltman coefficient 
functions 


In this section, the Feynman parameter integral repre¬ 
sentation of the Passarino-Veltman coefficient one, two, 
and three point functions are given. They are obtained 


Quantity 


Convention 


Metric signature g = diag(+, —, —, —) 

Spacetime dimension d = 4 — 2e 
Dirac matrix commutator = §[ 7 ^, 7 ,,] 

Fifth gamma matrix 75 = « 7 ° 7 1 7 2 7 3 
Chiral projectors Pl = §(1 — 75 ), Pr = §(1 + 75 ) 

Levi-civita symbol 


e 0123 = +1 


TABLE II. Conventions for spacetime quantities 
Function Diagram 


A 0 (m 0 ) 


m 0 


m 0 


B 0 (p 2 ,m 0 ,mi), 
and b^(p 2 ,mo,0) 


Co(pl,phq 2 ,m2,mi,m 0 ) , 
q 2 = {P 2 -P 1) 2 


m 1 


Pi 



TABLE III. Conventions for the arguments of the Passarino- 
Veltman functions 


[251 by writing tensor integrals as derivatives of the in¬ 
tegral representation of the corresponding scalar integral 
with respect to external momenta, and then matching the 
result to the respective covariant tensor decomposition. 

Ao^imo) = (4V) e ^^r(-l+e-r)mJ- e+r (Bl) 

2 r 


Bo...o„i...i.(p 2 ;mo,mi) = (4 irp 2 ) e 
dxx n 


j 2 +r+n 


/0 (p 2 x 2 + (—p 2 + m\ — itiq)x + Wq — ie) 


P(e-r) 
— (B2) 


( _ 1 \3+r+n 

!... 1 (P 2 - m ) = - 2r -P (1 + e — r) 


dxx n ( 1 — x) 


>0 (p 2 x 2 + (— p 2 + m 2 )x — ie) 


. \ l+e—r 


(B3) 
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c 0...0 1...1 2...2 (pi,pi,q 2 ;m 2 ,m 1 ,m 0 ) = ( 4 - 7 r/x 2 ) e 


^ 2) 3 + r + ni + n2 


r(l + e - r) 


x [ dy [ dzy ni z n2 [p\y 2 +p\z 2 + (-q 2 +p\+pl)yz+(-p\+ml-rnl)y + (-pl+ml-rnl)z + ml--i£\ 1 ' + ' 

Jo Jo 

(B4) 

The coefficient C function exhibits an invariance under the simultaneous interchange of indices n\ ■(->■ n 2 , external 
momenta p\ p 2 and internal masses m\ O m 2 , 


C 0...0 1...1 2...2 (pi,P 2 ,q 2 ;rn 2 ,mi,m 0 ) = C 0 ...o 1...1 2 ... 2 (p 2 ,pl,q 2 ',m 1 ,m 2 ,m 0 ) 


(B5) 


and is frequently employed during its reduction in the most general kinematic case (det Z ^ 0). 

In certain reduction formulae of C-functions, some terms are multiplied by e, which in the e —> 0 limit, pick up 
the UV-divergent parts of the coefficient functions in those terms. The UV divergent parts are readily obtained from 
the integral representation. They are controlled by the leading gamma function which for large enough r develops a 
1/e pole as e —> 0. When r is large, the integrand becomes polynomial in the Feynman parameters and are readily 
integrated with the help of the multinomial theorem. The needed UV-divergent parts are those of the B and C 
functions, shown below. 


B 0...0 1...1 (p 2 ;mo,mi) 


(- 1 )” 


E 


uv- 2 r B 

Div - ‘ k!+k 2 +k 3 =r 


a kl b k2 c k3 


1 


fci, k 2l k 3 J 2ki + k 2 + n + 1 e ’ 


where a = p , b = —p 2 + m\ — 7715 , and c = ?71 q, are polynomial coefficients of the integrand in (B2). 


(B 6 ) 


C 0...0 1...1 2 ... 2 {pl,pt,q 2 -,m 2 ,m 1 ,m 0 ) 

2r- ni n 2 

( — 1 ) ni + n 2 


UV- 

Div. 


2 r (r- 1)! 


E 


/ci + ...+fc6=r—1 


r ~ 1 ) n k ' h k2 r k \1 k *r k * f k ° ^ + fc 3 + fc 4 + rn)\(2k 2 + fc 3 + k 5 + n 2 )\ 1 . . 

ki,..., k 6 ) ( 2 ki + 2 k 2 + 2 k 3 + fc 4 + k 5 + ni + n 2 + 2 )! e ’ 


where a, b , c, d, e, and / are polynomial coefficients of the integrand in (B4) in the order displayed. 


Appendix C: Derivation of reduction formulae for C functions Case 2 


The derivation of the first equation in (191 begins with the Feynman parameter representation of the coefficient C 
function, 


C 0...0 1...1 2...2 (mo,s,m|;m 2 , 0 ,m o ) = ( 47 r/r 2 )' 


^_^^3+r-+ni+n 2 


r(l + e — r) 


f 1 f 1 ~ y _i_ , 

x dy dz y ni z n2 [m^y 2 + sz 2 + {—m 2 + TOq + s)yz + (—mg + m 2 — s)z + ml — *e] e r . (Cl) 

Jo Jo 


/ 0 ^0 

Upon making a change of integration variables y = 1 — y' and z = y'z' , the nested integrals are factored: 

7*1 p 1 

integrals = / dj/' 7 / ,_1+Il2 ~ 2£+2r (l — y ')" 1 / dV 2: ,n2 [s / 2 + (—s + m| — TOg)z 7 + mg — ie] 1 e+r . 
do Jo 


(C2) 


The y' integral gives the Euler Beta function, while the z' integral is identified as the integral representation of 
coefficient 13-function (B2). 


(-if 1 

2 


C 0...0 1...1 2...2 (mo,s,m|;m 2 ,0,m 0 ) = —B(n 2 - 2e + 2r,ni + 1)J3 0 ... 0 1...1 (s; m 0 , m 2 ) (C3) 
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As long as one of n 2 or r is non-zero, the Beta function is finite, and its expansion to 0(e) may be inserted yielding 
the first equation in (191. 

On the other hand, if n 2 = r = 0, the Beta function devel ops a 1/e pole, so that to 0(e), 

B(— 2e, ni + 1) = — H ni — e(-ff^ — Hn}). In this case, (Cl) is written as 


V -1-e 


1 f 1 

Ci.,.i (toq, s,ra|;m 2 ,0,m o ) = (47T+) e (-l) ni r(l + e)— / dz'(sz' 2 + z'(—s + m\ - ml) + ml - ie) 

ze Jo 

71 1 

+ ( 47 r/x 2 ) e (— l)" i r(l + e)(H ni + e(B[ 2 1 — H^)) [ dz'(sz’ 2 + z'{-s + m%-ml) + ml - ie)~ l e . (C4) 

Jo 

While the z' integral in the second line can be identified with the integral representation of the coefficient 13 function, 
the first line is identifiecf] as the integral representation of the scalar function Uo(mg, s, m 2 , 0 , mo) classified by 
Ellis and Zanderighi EH as IR-divergent triangle 6 . These identifications lead to the second equation of ( [19] ). 

If the off-shell momentum s is in the third argument, the derivation starts with the change of variables z = 1 — y — x 
in (B4) followed by an interchange of the x and y integrals to give 


0O...O„l...l„2...2( m 2) m O, s ; m O> m 2 ,O) = ( 47r M 2 )' 


^_l)3+r+ni+n2 


T(1 + e — r) 


/»1 /*1 — X 

x dx dy y ni (1 — x — y) n2 [mix 2 + sy 2 + (s + m§- ml)xy — 2m^x + (— s — ml + m%)y + ml — ie] 

Jo Jo 

The nested integrals are factored by making a further change of variables x = 1 — x' and y = y'x' to give 

/»1 rl 

integrals = / dx x /n i+ n 2 + 2 r-i- 2 e / ^ ^ [sy 2 + (—s + m\ — mfyy' + ttiq — ie] 1 e+r . 

J o Jo 


— 1 —e+r 


(C5) 


(C 6 ) 


The x 1 integral is straightforward. The y' integral can be brought to a recognizable form after expanding the factor 
(1 — y') n2 as a binomial series 

= ——-Tb-TEUlHn dy'y' ni+k [sy' 2 + (-s + m 2 2 -ml)y'+ ml-ie]~ l ~ e+r . (C7) 

ni + n 2 + 2r - 2e \ k J J 0 

The y' integral is now identified as the integral representation of the coefficient B function, yielding (20). 


Appendix D: Derivation of reduction formulae for C functions Case 4 


Two cases are distinguished for Case 4 (det Z = 0, Xq j = 0) depending on whether is vanishing. Although 
the steps below leading to (22) and (23) appear complicated, they essentially follow that of [12] for the evaluation of 


the scalar function 0 O . Beginning with the integral representation ( |B4[ ), a change of integration variables y = 1 — y' 
brings the Feynman integrals to the form 


C rv' 

integrals = / dy 1 dz (1 — y') ni z n2 [ay 12 + b z 2 + cy'z + dy'+ e z + f] 
Jo Jo 


— 1 —e+r 


(Dl) 


where a = p 2 , b = p|, c = q 2 — p\ — pi, d = —p 2 + ml — mf, e = pf — q 2 — ml + m 2 , and / = m 2 — ie. 

Under the assumption that p\^ 0, a second change of variables is made z = z' + ay' , with a = + chosen to make 
the coefficient of y' 2 in square brackets vanish. 


integrals = f dy' f dz' (1 — y') ni (z' + ay') n2 [bz' 2 + {c + 2ba)y'z' + (d+ ea)y' + ez' + f] 

Jo J—av 


— 1— e+r 


(D2) 


4 see http://qcdloop.fnal.gov/tridiv6.pdf 
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The choice for a implies that c + 2 ba vanishes, and the kinematic relations det Z = Xqj = 0 imply that d + ea 
vanishes, yielding 

nl n(l — Ot)y 

integrals = / dy dz (1 - y) n '{z + ay) n *[bz 2 + e * + /] _1_e+r , (D3) 

J 0 J— ay 

where the primes have been omitted. The binomial theorem is applied to the factor ( z+ay ) n2 = Yj (” 2 ) a n ' 2 ~ 3 y n<2 ~ 3 z 3 , 
and the order of integrations is interchanged so that 


itegrals = ^ 
l=o 


r-l-a 


dz 


L J 0 


/z/(l-a) 


dy - 


dz 


' —z j a 


dy (1 — y) ni y n2 3 z 3 [bz 2 + ez + f] 


• (D4) 


The y integrals in both terms yield terminating hypergeometric series most compactly written in terms of the incom¬ 
plete Beta function: 



dy( 1 - y) ni y n2 - j 


ni\(n 2 - j)! 

(ni + n 2 - j + 1 )! 

ni!(n 2 - j)! 

(ni + n 2 - j + 1 )! 


B x (n 2 -j + l,ni + 1 ) 


(_^kj£n 2 —j+k+l 

^ (n 2 -j + k + 1 ) 



(D5) 


Since the first term of (D51 is common to both integrations in (D4), they are combined to yield a total of three terms: 


n 2 


integrals = / 
l=o 




ni\(n 2 - j)\ f 1 

L(m +n 2 -j + 1)! J_ c 


dz 


z n 2 +i 


(bz 


f) 


1+e —r 


nl ° +/(!-<*)+2 + Urn + f)^ 


(bz 


f) 


1+e—r 


dz 


B - z/a (n 2 -j + l,m + 1 ) z 3 


(bz 


f) 


1+e —r 


(D 6 ) 


A change of integration variables is carried out in each term to stretch their ranges to 0 —1 1: In the first integral, 
z = z' — a, in the second integral z = (1 — a)z ', and in the third integral 2 = —az'. Consequently, the polynomials 
bz 2 + ez + f take the shape of integrands for the B function^ 


First term: 
Second term: 
Third term: 


p 2 z' 2 + (~P 2 + m 2 — rrifyz' + rriQ — is 
q 2 z’ 2 + (—q 2 + m 2 — m 2 )z' + m\ — is 
p\z' 2 + (—p\ + TOp — m 2 )z' + m\—is 


= PzW 
= Pl2(z') 
= Pi(z') 


Upon inserting the series representation of the incomplete Beta function (D5) the result is (after dropping the primes 
on z ) 


integrals = 

l=o 


(- + 

2—1 r, 0 _ n 




n 1 


rci!(n 2 - j)! 

(ni +n 2 - j + 1 )! 


dz(z — a) 3 P 2 (z) 


— 1 — e+r 


k —0 


n 2 —j + k + l\ k 


- (1 - a) 3+1 / dzz n2+k+1 P 12 (z)- l ~ t+r + (~a) 3+1 / d 2 2 " 2 +?c+ 1 P 1 ( 2 )- 1 - e+r 


(D7) 


In the first term, the binomial theorem is applied to (z — a) 3 = Yk (1) (~ a ) 3 ~ k z k , and the three 2 integrals are finally 
identified as integral representations of the coefficient B functions upon which (22) is obtained. 

If p 2 = 0, equation (22) breaks down and another formulae is needed. In this c ase, det Z = 0 implies p\ = q 2 and 
X 0 j = 0 implies m 0 = m 2 provided pf ^ 0. With these relations, the integrals in (D1) are already factored: 


r 1 rv‘ 

integrals = / dy 1 dz (1 — y') ni z n2 [pfy 12 + (— p 2 + mg — m 2 )y' + m\ — is\ 

Jo Jo 


— 1—e+r 


(D 8 ) 


5 These quadratic polynomials may be recognized as the ‘pinch 
functions’ originating from the three cut channels of the triangle 


graph. 
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The 2 integration gives a factor l/(ri 2 +1), and the factor (1 — y') ni = J2k (T) ( — y) k is expressed as a binomial series. 


1 1 / \ r 1 

integrals = n ^ + l X] t -1 )* J d y'y ,n2+k+1 [ply 12 + (- p\ + m 2 0 - ml)y’ + m\- ie\ 


— 1 — e+r* 


(D9) 


Eqn (23) is obtained after identifying the y' integration as the integral representation of the coefficient B function. If 
all external invariants are vanishing p\ = p% = q 2 = 0 then neither (22) nor (23) are valid, and Case 5 is needed. 


Appendix E: Coefficient B functions with r = — 1 


The two new reduction algorithms for C functions 
(Cases 2 and 4) require extending the set of basis func¬ 
tions to include B functions in which the index r in (B2) 
continued to — 1 . 


is continued to — 1 . In a certain sense, these new re¬ 
duction formulae may be closely related to those in [26] . 
There, the authors present different reduction formulae 
for coefficient functions which likewise require extending 
the set of basis functions to scalar functions with repeated 
propagators. 

A set of explicit expressions for the general case and at 
singular points is constructed and included in the Pack¬ 
age -X source file (see Table]]]). The integration is straight¬ 
forward in most cases. The functions are UV-finite for 
all n > 0 , but with many kinematic configurations devel¬ 
oping IR-divergent 1/e poles. 

However, there are three kinematic cases, all corre¬ 
sponding to physical threshold for which the Feynman 
parameter integral nominally diverges even for finite but 
infinitesimal e. To handle these cases, e is taken suffi¬ 
ciently large and negative so that the integral converges, 
and then analytically continued to e —> 0. The results of 
these integrations are given below: 


B o...o 1...1 (to 2 ; 0 , mi) 


- 2 , 47 T/r 2 ,e f 1 

= r (! + e ) / dx ‘ 

,n 1 Jo 


2—2e 


(El) 


r 


B o...o 1...1 (to o! toq, 0) 


= ~2 (~^~) e ( — 1)" +1 T(1 + e) C dxx n {x- l)- 2 - 
TOq Toq J 0 

+ln(^r) + 2 F n _i - 2 ^, n > 1 

n = 0 


(E2) 


-B 0...0 1...1 ((to.o+toi) 2 ;to 0 ,toi) = 


2 (— 1 ) 


n+1 


(mo+TOi) 2 


47T/i 2 


0 r(1+e) i 


dx- 


(too+toi) 2 ^ 1 '’Jo ^ [(cc — rr + ) 2 ] 1+e 

[ Y°b— +„*r 1 ln (1-^L) 

L k + 1 + 1 x + ’ 

too 


2(_l) rl + 1 r2_ 2 T n ~ 

(m 0 +m^) 


k =0 


1 


y n,0 


1 — x + x + 


X+ = 


mo — mi 


(E3) 


Among these integrals, only ( |E3[ ) gives numerical results 
that are related to limiting values as threshold is reached: 
the real part of (E3) corresponds to the limiting value of 


ReH(s; too, toi) when approached from above threshold, 
and the imaginary part corresponds to the limiting value 
of ImI3(s; too, toi) when approached below threshold. 

That the integrals ( E1||E3 ) give numerical results that 
do not match their limiting values as threshold is reached 
are not of any concern. The results above should be 
viewed as ill-defined divergent integrals arising at inter¬ 
mediate stages in the reduction of coefficient C functions. 
They serve to facilitate the analytic cancellation of these 
integrals at the end of a physically meaningful computa¬ 
tion, such as for the electromagnetic contribution to the 
electron anomalous magnetic moment. 
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